#include <iostream>
#include <pthread.h>
#include <malloc.h>
#include <time.h>
#include <stdlib.h>
#include <semaphore.h>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX¡¢AVX2

using namespace std;
const int NUM_THREADS=5;
int n=1000;
float **a;
float **b1;
float **b2;
float **c1;
float **c2;
float **c3;
float **d;
float **e;
void starta(float** x, int n) {//初始化数组
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < i; j++)
            x[i][j] = 0;
        x[i][i] = 1.0;
        for (int j = i + 1; j < n; j++)
            x[i][j] = rand();
    }
    for (int k = 0; k < n; k++) {
        for (int i = k + 1; i < n; i++) {
            for (int j = 0; j < n; j++) {
                x[i][j] += x[k][j];
            }
        }
    }
}
void clone(float **A, float** B,int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = A[i][j];
        }
    }
}
void gaussionone(float** x, int n) {//高斯消元串行算法
    for (int k = 0; k < n; k++) {
        for (int j = k + 1; j < n; j++) {
            x[k][j] /= x[k][k];
        }
        x[k][k] = 1.0;
        for (int i = k + 1; i < n; i++) {
            for (int j = k + 1; j < n; j++) {
                x[i][j] -= x[i][k] * x[k][j];
            }
            x[i][k] = 0;
        }
    }
}
//动态线程版本
typedef struct{
    int k;//消去的轮数
    int t_id;//线程id
}threadParam_t;
void *threadFunc(void *param){
    threadParam_t *p=(threadParam_t*)param;
    int k=p->k;
    int t_id=p->t_id;
    int i=k+t_id+1;
    for(int j=k+1;j<n;j++){
        b1[i][j]-=b1[i][k]*b1[k][j];
    }
    b1[i][k]=0;
    pthread_exit(NULL);
}
void* threadFuncsse(void* param)
{
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;
    int t_id = p->t_id;
    int i = k + t_id + 1;
    __m128 vaik = _mm_set1_ps(b2[i][k]);
    int j;
    for (j = k + 1; j + 4 <= n; j += 4)
    {
        __m128 vakj = _mm_loadu_ps(&b2[k][j]);
        __m128 vaij = _mm_loadu_ps(&b2[i][j]);
        __m128 vx = _mm_mul_ps(vakj, vaik);
        vaij = _mm_sub_ps(vaij, vx);
        _mm_storeu_ps(&b2[i][j], vaij);
    }
    for (; j < n; j++)
    {
        b2[i][j] = b2[i][j] - b2[k][j] * b2[i][k];
    }
    b2[i][k] = 0;
    pthread_exit(NULL);
}
void gaussionaction(float** x,int n){
    for(int k=0;k<n;k++){
        for(int j=k+1;j<n;j++){
            x[k][j]/=x[k][k];
        }
        x[k][k]=1.0;
        int worker_count=n-1-k;
        pthread_t* handles = (pthread_t*)malloc(worker_count * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(worker_count * sizeof(threadParam_t));
        for(int t_id=0;t_id<worker_count;t_id++){
            param[t_id].k=k;
            param[t_id].t_id=t_id;
        }
        for(int t_id=0;t_id<worker_count;t_id++){
            pthread_create(&handles[t_id], NULL, threadFunc, (void*)&param[t_id]);
        }
        for(int t_id=0;t_id<worker_count;t_id++){
            pthread_join(handles[t_id], NULL);
        }
    }
}
void gaussionactionsse(float** x,int n)
{
    for (int k = 0; k < n; k++)
    {
        __m128 vt = _mm_set1_ps(x[k][k]);
        int j;
        for (j = k + 1; j + 4 <= n; j += 4)
        {
            __m128 va = _mm_loadu_ps(&x[k][j]);
            va = _mm_div_ps(va, vt);
            _mm_storeu_ps(&x[k][j], va);
        }
        for (; j < n; j++)
        {
            x[k][j] = x[k][j] / x[k][k];
        }
        x[k][k] = 1.0;

        int worker_count = n - k - 1;
        pthread_t* handles = (pthread_t*)malloc(worker_count * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(worker_count * sizeof(threadParam_t));
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            pthread_create(&handles[t_id], NULL, threadFuncsse, (void*)&param[t_id]);
        }
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            pthread_join(handles[t_id], NULL);
        }
    }
}

//静态线程+信号量版本
typedef struct{
    int t_id;//线程id
}threadParam_t1;
sem_t sem_main;
sem_t sem_workerstart[NUM_THREADS];
sem_t sem_workerend[NUM_THREADS];
void *threadFunc1(void *param){
    threadParam_t1 *p=(threadParam_t1*)param;
    int t_id=p->t_id;
    for(int k=0;k<n;k++){
        sem_wait(&sem_workerstart[t_id]);
        for(int i=k+1+t_id;i<n;i+=NUM_THREADS){
            for(int j=k+1;j<n;j++){
                c1[i][j]-=c1[i][k]*c1[k][j];
            }
            c1[i][k]=0.0;
        }
        sem_post(&sem_main);
        sem_wait(&sem_workerend[t_id]);
    }
    pthread_exit(NULL);
}
void gaussionquite(float** x,int n){
    sem_init(&sem_main,0,0);
    for(int i=0;i<NUM_THREADS;i++){
        sem_init(&sem_workerstart[i],0,0);
        sem_init(&sem_workerend[i],0,0);
    }
    pthread_t handles[NUM_THREADS];
    threadParam_t1 param[NUM_THREADS];
    for(int t_id=0;t_id<NUM_THREADS;t_id++){
            param[t_id].t_id=t_id;
            pthread_create(&handles[t_id], NULL, threadFunc1, (void*)&param[t_id]);
    }

    for(int k=0;k<n;k++){
        for(int j=k+1;j<n;j++){
            x[k][j]/=x[k][k];
        }
        x[k][k]=1.0;
        for(int t_id=0;t_id<NUM_THREADS;t_id++){
            sem_post(&sem_workerstart[t_id]);
        }
        for(int t_id=0;t_id<NUM_THREADS;t_id++){
            sem_wait(&sem_main);
        }
        for(int t_id=0;t_id<NUM_THREADS;t_id++){
            sem_post(&sem_workerend[t_id]);
        }

    }
    for(int t_id=0;t_id<NUM_THREADS;t_id++){
        pthread_join(handles[t_id], NULL);
    }
    sem_destroy(sem_workerstart);
    sem_destroy(sem_workerend);
}
void* threadFunc1sse(void* param)
{
    threadParam_t1* p = (threadParam_t1*)param;
    int t_id = p->t_id;
    for (int k = 0; k < n; k++)
    {
        sem_wait(&sem_workerstart[t_id]);
        for (int i = k + t_id + 1; i < n; i +=NUM_THREADS)
        {
            __m128 vaik = _mm_set1_ps(c2[i][k]);
            int j;
            for (j = k + 1; j + 4 <= n; j += 4)
            {
                __m128 vakj = _mm_loadu_ps(&c2[k][j]);
                __m128 vaij = _mm_loadu_ps(&c2[i][j]);
                __m128 vx = _mm_mul_ps(vakj, vaik);
                vaij = _mm_sub_ps(vaij, vx);
                _mm_storeu_ps(&c2[i][j], vaij);
            }
            for (; j < n; j++)
            {
                c2[i][j]-= c2[k][j] * c2[i][k];
            }
            c2[i][k] = 0;
        }
        sem_post(&sem_main);
        sem_wait(&sem_workerend[t_id]);
    }
    pthread_exit(NULL);
}
void gaussionquitesse(float** x,int n)
{
    sem_init(&sem_main, 0, 0);
    for (int i = 0; i <NUM_THREADS; i++)
    {
        sem_init(&sem_workerstart[i], 0, 0);
        sem_init(&sem_workerend[i], 0, 0);
    }
    pthread_t handles[NUM_THREADS];
    threadParam_t1 param[NUM_THREADS];
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        param[t_id].t_id = t_id;
        pthread_create(&handles[t_id], NULL, threadFunc1sse, (void*)&param[t_id]);
    }
    for (int k = 0; k < n; k++)
    {
        __m128 vt = _mm_set1_ps(x[k][k]);
        int j;
        for (j = k + 1; j + 4 <= n; j += 4)
        {
            __m128 va = _mm_loadu_ps(&x[k][j]);
            va = _mm_div_ps(va, vt);
            _mm_storeu_ps(&x[k][j], va);
        }
        for (; j < n; j++)
        {
            x[k][j] /= x[k][k];
        }
        x[k][k] = 1.0;
        for (int t_id = 0; t_id <NUM_THREADS; t_id++)
        {
            sem_post(&sem_workerstart[t_id]);
        }
        for (int t_id = 0; t_id <NUM_THREADS; t_id++)
        {
            sem_wait(&sem_main);
        }
        for (int t_id = 0; t_id <NUM_THREADS; t_id++)
        {
            sem_post(&sem_workerend[t_id]);
        }
    }
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        pthread_join(handles[t_id], NULL);
    }
    sem_destroy(&sem_main);
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        sem_destroy(&sem_workerstart[t_id]);
        sem_destroy(&sem_workerend[t_id]);
    }
}
void* threadFunc1avx(void* param)
{
    threadParam_t1* p = (threadParam_t1*)param;
    int t_id = p->t_id;
    for (int k = 0; k < n; k++)
    {
        sem_wait(&sem_workerstart[t_id]);
        for (int i = k + t_id + 1; i < n; i +=NUM_THREADS)
        {
            __m256 vaik = _mm256_set1_ps(c3[i][k]);
            int j;
            for (j = k + 1; j + 8 <= n; j += 8)
            {
                __m256 vakj = _mm256_loadu_ps(&c3[k][j]);
                __m256 vaij = _mm256_loadu_ps(&c3[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_storeu_ps(&c3[i][j], vaij);
            }
            for (; j < n; j++)
            {
                c3[i][j]-= c3[k][j] * c3[i][k];
            }
            c3[i][k] = 0;
        }
        sem_post(&sem_main);
        sem_wait(&sem_workerend[t_id]);
    }
    pthread_exit(NULL);
}
void gaussionquiteavx(float** x,int n)
{
    sem_init(&sem_main, 0, 0);
    for (int i = 0; i <NUM_THREADS; i++)
    {
        sem_init(&sem_workerstart[i], 0, 0);
        sem_init(&sem_workerend[i], 0, 0);
    }
    pthread_t handles[NUM_THREADS];
    threadParam_t1 param[NUM_THREADS];
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        param[t_id].t_id = t_id;
        pthread_create(&handles[t_id], NULL, threadFunc1avx, (void*)&param[t_id]);
    }
    for (int k = 0; k < n; k++)
    {
        __m256 vt = _mm256_set1_ps(x[k][k]);
        int j;
        for (j = k + 1; j + 8 <= n; j += 8)
        {
            __m256 va = _mm256_loadu_ps(&x[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_storeu_ps(&x[k][j], va);
        }
        for (; j < n; j++)
        {
            x[k][j] /= x[k][k];
        }
        x[k][k] = 1.0;
        for (int t_id = 0; t_id <NUM_THREADS; t_id++)
        {
            sem_post(&sem_workerstart[t_id]);
        }
        for (int t_id = 0; t_id <NUM_THREADS; t_id++)
        {
            sem_wait(&sem_main);
        }
        for (int t_id = 0; t_id <NUM_THREADS; t_id++)
        {
            sem_post(&sem_workerend[t_id]);
        }
    }
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        pthread_join(handles[t_id], NULL);
    }
    sem_destroy(&sem_main);
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        sem_destroy(&sem_workerstart[t_id]);
        sem_destroy(&sem_workerend[t_id]);
    }
}
//静态线程+信号量+三重循环纳入线程函数
typedef struct{
    int t_id;
}threadParam_t2;
sem_t sem_leader;
sem_t sem_Divsion[NUM_THREADS-1];
sem_t sem_Elimination[NUM_THREADS-1];
void *threadFunc2(void *param){
    threadParam_t2 *p=(threadParam_t2*)param;
    int t_id=p->t_id;
    for (int k = 0; k < n; k++)
    {
        if (t_id == 0)
        {
            for (int j = k + 1; j < n; j++)
            {
                d[k][j] = d[k][j] / d[k][k];
            }
            d[k][k] = 1.0;
        }
        else
        {
            sem_wait(&sem_Divsion[t_id - 1]);
        }
        if (t_id == 0)
        {
            for (int i = 0; i <NUM_THREADS- 1; i++)
            {
                sem_post(&sem_Divsion[i]);
            }
        }
        for (int i = k + t_id + 1; i < n; i +=NUM_THREADS)
        {
            for (int j = k + 1; j < n; j++)
            {
                d[i][j] = d[i][j] - d[i][k] * d[k][j];
            }
            d[i][k] = 0.0;
        }
        if (t_id == 0)
        {
            for (int i = 0; i <NUM_THREADS- 1; i++)
            {
                sem_wait(&sem_leader);
            }
            for (int i = 0; i <NUM_THREADS- 1; i++)
            {
                sem_post(&sem_Elimination[i]);
            }
        }
        else
        {
            sem_post(&sem_leader);
            sem_wait(&sem_Elimination[t_id - 1]);
        }
    }
    pthread_exit(NULL);
}
void gaussionquitesan()
{
    sem_init(&sem_leader, 0, 0);
    for (int i = 0; i < NUM_THREADS- 1; i++)
    {
        sem_init(&sem_Divsion[i], 0, 0);
        sem_init(&sem_Elimination[i], 0, 0);
    }
    pthread_t handles[NUM_THREADS];
    threadParam_t2 param[NUM_THREADS];
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        param[t_id].t_id = t_id;
        pthread_create(&handles[t_id], NULL, threadFunc2, (void*)&param[t_id]);
    }
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        pthread_join(handles[t_id], NULL);
    }
    sem_destroy(&sem_leader);
    for (int t_id = 0; t_id <NUM_THREADS - 1; t_id++)
    {
        sem_destroy(&sem_Divsion[t_id]);
        sem_destroy(&sem_Elimination[t_id]);
    }
}
//静态线程+barrier
typedef struct{
    int t_id;
}threadParam_t3;
pthread_barrier_t barrier_Divsion;
pthread_barrier_t barrier_Elimination;
void *threadFunc3(void *param){
    threadParam_t3 *p=(threadParam_t3*)param;
    int t_id=p->t_id;
    for(int k=0;k<n;k++){
        if(t_id==0){
            for(int j=k+1;j<n;j++){
                e[k][j]/=e[k][k];
            }
            e[k][k]=0;
        }
        pthread_barrier_wait(&barrier_Divsion);
        for(int i=k+1+t_id;i<n;i+=NUM_THREADS){
            for(int j=k+1;j<n;j++){
                e[i][j]-=e[i][k]*e[k][j];
            }
            e[i][k]=0.0;
        }
        pthread_barrier_wait(&barrier_Elimination);
    }
    pthread_exit(NULL);
}
void gaussionquitebarrier(){
    pthread_barrier_init(&barrier_Divsion, NULL,NUM_THREADS);
    pthread_barrier_init(&barrier_Elimination, NULL,NUM_THREADS);

    pthread_t handles[NUM_THREADS];
    threadParam_t3 param[NUM_THREADS];
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        param[t_id].t_id = t_id;
        pthread_create(&handles[t_id], NULL, threadFunc3, (void*)&param[t_id]);
    }
    for (int t_id = 0; t_id <NUM_THREADS; t_id++)
    {
        pthread_join(handles[t_id], NULL);
    }
    pthread_barrier_destroy(&barrier_Divsion);
    pthread_barrier_destroy(&barrier_Elimination);
}
int main()
{
   for (int k = 100; k <2000; k += 100) {
    n = k;
    a=new float*[n];
    b1=new float*[n];
    b2=new float*[n];
    c1=new float*[n];
    c2=new float*[n];
    c3=new float*[n];
    d=new float*[n];
    e=new float*[n];
    for(int i=0;i<n;i++){
        a[i]=new float[n];
        b1[i]=new float[n];
        b2[i]=new float[n];
        c1[i]=new float[n];
        c2[i]=new float[n];
        c3[i]=new float[n];
        d[i]=new float[n];
        e[i]=new float[n];
    }
    starta(a,n);
    clone(a,b1,n);
    clone(a,b2,n);
    clone(a,c1,n);
    clone(a,c2,n);
    //clone(a,c3,n);
    clone(a,d,n);
    clone(a,e,n);

    clock_t start1,finish1;
    start1 = clock();
    gaussionone(a, n);
    finish1 = clock();
    cout <<n<<"串行"<<(double)(finish1 - start1) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start2,finish2;
    start2 = clock();
    gaussionaction(b1, n);
    finish2 = clock();
    cout <<n<<"动态线程"<<(double)(finish2 - start2) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start3,finish3;
    start3 = clock();
    gaussionactionsse(b2, n);
    finish3 = clock();
    cout <<n<<"动态线程+sse"<<(double)(finish3 - start3) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start4,finish4;
    start4 = clock();
    gaussionquite(c1, n);
    finish4 = clock();
    cout <<n<<"静态线程+信号量"<<(double)(finish4 - start4) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start5,finish5;
    start5 = clock();
    gaussionquitesse(c2, n);
    finish5 = clock();
    cout <<n<<"静态线程+信号量+sse"<<(double)(finish5 - start5) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start6,finish6;
    start6 = clock();
    gaussionquiteavx(c3, n);
    finish6 = clock();
    cout <<n<<"静态线程+信号量+avx"<<(double)(finish6 - start6) * 1000 / CLOCKS_PER_SEC << "ms" << endl;
    /*clock_t start7,finish7;
    start7 = clock();
    gaussionquitesse(c2, n);
    finish7 = clock();
    cout <<n<<"静态线程+信号量+col"<<(double)(finish7 - start7) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start8,finish8;
    start8 = clock();
    gaussionquitesse(c2, n);
    finish8 = clock();
    cout <<n<<"静态线程+信号量+row"<<(double)(finish8 - start8) * 1000 / CLOCKS_PER_SEC << "ms" << endl;*/

    clock_t start9,finish9;
    start9 = clock();
    gaussionquitesan();
    finish9 = clock();
    cout <<n<<"静态线程+信号量+三重循环"<<(double)(finish9 - start9) * 1000 / CLOCKS_PER_SEC << "ms" << endl;

    clock_t start,finish;
    start = clock();
    gaussionquitebarrier();
    finish = clock();
    cout <<n<<"静态线程+barrier"<<(double)(finish - start) * 1000 / CLOCKS_PER_SEC << "ms" << endl;
    cout<<endl;
}
}
